Comparison Theorems for Stochastic Chemical Reaction Networks

Continuous-time Markov chains are frequently used as stochastic models for chemical reaction networks, especially in the growing field of systems biology. A fundamental problem for these Stochastic Chemical Reaction Networks (SCRNs) is to understand the dependence of the stochastic behavior of these systems on the chemical reaction rate parameters. Towards solving this problem, in this paper we develop theoretical tools called comparison theorems that provide stochastic ordering results for SCRNs. These theorems give sufficient conditions for monotonic dependence on parameters in these network models, which allow us to obtain, under suitable conditions, information about transient and steady-state behavior. These theorems exploit structural properties of SCRNs, beyond those of general continuous-time Markov chains. Furthermore, we derive two theorems to compare stationary distributions and mean first passage times for SCRNs with different parameter values, or with the same parameters and different initial conditions. These tools are developed for SCRNs taking values in a generic (finite or countably infinite) state space and can also be applied for non-mass-action kinetics models. When propensity functions are bounded, our method of proof gives an explicit method for coupling two comparable SCRNs, which can be used to simultaneously simulate their sample paths in a comparable manner. We illustrate our results with applications to models of enzymatic kinetics and epigenetic regulation by chromatin modifications. Supplementary Information The online version contains supplementary material available at 10.1007/s11538-023-01136-5.

Theorem S.1. Let X be an irreducible continuous-time Markov chain with state space X and infinitesimal generator Q. Suppose V : X → R + is norm-like, that is {x ∈ X : V (x) ≤ a} is compact 6 for each a ∈ R + . Further assume that for some c > 0, d > 0 and a compact set C, QV (x) ≤ −c + d1 C (x), for all x ∈ X . (S.1) Then, X is non-explosive and positive recurrent, and has a unique stationary distribution π. If instead of (S.1), we have that for some c ′ > 0 and d ′ > 0, for all x ∈ X , (S. 2) then (S.1) automatically holds and the consequences stated above hold, and in addition, the stationary distribution satisfies and there is 0 < B < ∞ and 0 < β < 1 such that for all t ≥ 0 and x ∈ X , where P xy (t) = P[X(t) = y|X(0) = x] and ∥P x• (t) − π∥ V +1 = sup |g|≤V +1 y∈X (P xy (t) − π y )g(y) .
Remark S.1. When the second inequality in (S.3) holds, we say that the Markov chain is exponentially ergodic in the (V + 1)-norm.
Proof. We will verify the sufficient conditions for each of non-explosion, positive recurrence and exponential ergodicity given in Meyn & Tweedie [1]. Note that X is a Borel rightprocess under the definition in Sharpe [2]. In addition, since the state space is discrete, each compact set is finite and therefore petite. For m ∈ Z + , if Q m is the infinitesimal generator for the Markov chain X killed upon exit from O m = {x ∈ X : V (x) ≤ m} 7 , then Q m V (x) ≤ QV (x) for x ∈ O m . It then follows from (S.1) that conditions (CD0) and (CD2) (with f = 1) in [1] hold with Q m in place of A m there. By Theorem 2.1 and Theorem 4.2 in [1], the Markov chain is non-explosive and positive recurrent, and it has a unique stationary distribution π.
On the other hand, if (S.2) holds, then (S.1) holds, using the norm like property of V . Furthermore, (S.2) implies that conditions (CD0), (CD2) (with f = V + 1) and (CD3) in [1] hold with Q m in place of A m there. By Theorem 2.1, Theorem 4.2 and Theorem 6.1 in [1], the Markov chain is non-explosive and positive recurrent, with a unique stationary distribution π such that π(V ) < ∞, and it is exponentially ergodic in the (V + 1)-norm, that is the second inequality in (S.3) holds for all t ≥ 0 and x ∈ X . For fixed t ≥ 0 and x ∈ X , setting g(y) = sgn(P xy (t) − π y ), for y ∈ X , we have that |g| ≤ 1 ≤ V + 1, and yielding the first inequality in (S.3).

S.1.2 Application to Example 4.2
For Example 4.2, we first show that the Markov chain is irreducible. For this, consider x • = (0, 0, E tot , 0) and any fixed state Starting at x • , by having reaction 5 ○ fire x 1 + E tot − x 3 + x 2 times in succession, then having reaction 1 ○, immediately followed by reaction 3 ○, fire x 2 times in succession and then reaction 1 ○ fire E tot − x 3 times, without any other reactions firing, we see that the Markov chain can transition with positive probability from x • to x. Since each reaction is reversible, it also follows that the Markov chain can transition from x to x • with positive probability. Thus, the Markov chain is irreducible. Next we will introduce a norm-like function V and show that (S.1) holds. For each x ∈ X , let Notice that b > 0 since E tot ≥ 1. Then, for each a ∈ R + , {x ∈ X : V (x) ≤ a} consists of finitely many states, and for any x ∈ X , The last equality uses the fact that x 3 = E tot − x 4 . For the following, we note that −2bx 2 3 + (2b − 1)x 3 = (−2bx 3 + (2b − 1))x 3 ≤ −1 when x 3 ≥ 1. We now consider two cases for QV (x): when x 3 = 0 and x 3 > 0. For the first case, when x 3 = 0, we have x 4 = E tot and As a quadratic function, the last expression is bounded above by since b is chosen as in (S.4). For the second case when x 3 ∈ {1, . . . , E tot }, we have where we have used the fact that 2ab ≤ a 2 + b 2 for the second inequality and the last expression will be less than or equal to −1 whenever x 1 or x 2 is sufficiently large.
Let C = {x ∈ X : QV (x) > −1}. Then, by the above, C consists of finitely many points, which implies that C is a compact set. Then, (S.1) holds with c = 1 and d = 1 + max x∈C QV (x) ∨ 0. It follows by Theorem S.1 that the Markov chain in Example 4.2 is non-explosive and positive recurrent, and has a unique stationary distribution π.

S.1.3 Application to Example 4.5
For Example 4.5, we first check that the Markov chain is irreducible. For this, consider x • = (0, D tot , 0) and any fixed state Starting at x • , by having reaction 5 ○ fire x 3 times, then reaction 3 ○ fire D tot − x 2 times, and finally, reaction 2 ○ fire x 1 times, without any other reactions firing, we see that the Markov chain can transition from x • to x with positive probability. For the reverse transition, by having reaction 6 ○ fire x 3 times, then reaction 4 ○ fire x 1 times, and finally, reaction 1 ○ fire D tot −x 2 times, we see that the Markov chain can transition from x to x • with positive probability. Thus, the Markov chain is irreducible. Next we will introduce a norm-like function V and show that (S.2) holds. For each x 3 ≤ a} consists of finitely many states for each a ∈ R + , and for any x ∈ X , where c ′ = κ 6a and d ′ = κ 5a . Therefore, we conclude by Theorem S.1 that the Markov chain in Example 4.5 is non-explosive and positive recurrent with a unique stationary distribution π such that π(V ) < ∞, and it is exponentially ergodic in the (V + 1)-norm.

S.2.1 Example 4.2
The set of reactions associated to the chemical reaction system in Fig. 3

S.2.2 Example 4.3
The set of five reactions associated to the chemical reaction system in Fig. 4(a) is given by

S.3 A generalization of Theorem 3.3
Here, we provide a more general version of Theorem 3.3. The simpler form given as Theorem 3.3 in the main text was used there because it is more straightforward to state and the conditions are easier to verify. However, the more general version of the theorem provided in this section can be useful in some cases, such as Example 4.3 (see Section S.3.2). The generalization relies on the idea that grouping of vectors can be more general than what is described in (3.12), and so we introduce the following assumption.
Consider a non-empty set X ⊆ Z d + , suppose that Assumption S.1 holds and consider two collections of non-negative functions on X , Υ = (Υ 1 , . . . , Υ n ) andΥ = (Υ 1 , . . . ,Υ n ), such that (3.3) holds and the associated continuous-time Markov chains do not explode in finite time. Further suppose that both of the following conditions hold: Then, for each pair x • ,x • ∈ X such that x • ≼ Ax • , there exists a probability space (Ω, F, P) with two continuous-time Markov chains X = {X(t), t ≥ 0} andX = {X(t), t ≥ 0} defined there, each having state space X ⊆ Z d + , with infinitesimal generators Q andQ, associated with Υ andΥ respectively, with initial conditions X(0) = x • andX(0) =x • and such that: Furthermore, in either case, for q * < q ≤ p k , we have ⟨A i• , v σ(q) ⟩ = 0.
The proof for Theorem S.2 can be found in Section S.3.1. With Theorem S.2 in place, we can extend Theorems 3.4 and 3.5 by adding an additional alternative condition (iv): Assumption S.1 holds, and conditions (i) and (ii) in Theorem S.2 are satisfied. If this is satisfied instead of one of (i) -(iii) in Theorems 3.4 and 3.5, then the conclusions of these theorems about (mean) first passage times and stationary distributions will still hold.

S.3.1 Proof of Theorem S.2
We initially assume that sup x∈X Υ j (x) < ∞ and sup x∈XΥj (x) < ∞ for every 1 ≤ j ≤ n, and let λ > 0 such that (5.2) holds. We shall relax these assumptions later. Further suppose that Assumption S.1 and condition (i) of Theorem S.2 both hold. For x ∈ X , define I k (x), I k q (x), Ψ λ ,Ȋ k (x),Ȋ k q (x) andΨ λ in the same manner as in the proof of Theorem 3.3 (see (5.32) -(5.34)), with {G k | 1 ≤ k ≤ s}, {p k | 0 ≤ k ≤ s} and σ as in Assumption S.1. Our proof of Theorem S.2 has some elements that are the same as those for the proof of Theorem 3.3. However, some additional elements are needed. We give the details for completeness. As for Theorem 3.3, Ψ λ (·, ·) andΨ λ (·, ·) are well-defined as X -valued functions.

(S.8)
Proof. First, we note that Ψ λ ,Ψ λ have the following properties: for every u ∈ We also have that, for 1 ≤ k ≤ s and j ∈ G k , x ≼ A y + v j if and only if To prove (S.8), we first consider the situation where y ∈ int(K A + x) = {w ∈ R d | Ax < Aw}. Then, for each 1 ≤ i ≤ m, ⟨A i• , y − x⟩ > 0 and since A ∈ Z m×d and y − x ∈ Z d , we have ⟨A i• , y − x⟩ ≥ 1. This implies that for each 1 ≤ k ≤ s and j ∈ G k , In addition, for 1 ≤ k ≤ s and j, ℓ ∈ G k , It follows from (S.16) -(S.18) that if y ∈ int(K A + x) ∩ X , then for any 1 ≤ k ≤ s and j, ℓ ∈ G k : x We also have, by assumption, that x ≼ A y. It follows that if y ∈ int(K A + x) ∩ X , then {x, x + v ℓ | ℓ ∈ G k } ≼ A {y, y + v j | j ∈ G k } for 1 ≤ k ≤ s and consequently (S.8) holds for all u ∈ [0, 1]. Now, we turn to the other situation where y ∈ ∂ i (K A + x) ∩ X for some 1 ≤ i ≤ m. Then K y := {i | ⟨A i• , y⟩ = ⟨A i• , x⟩, 1 ≤ i ≤ m} is non-empty. Let u ∈ [0, 1]. We consider two cases.
Fix such an index j. Consider the unique 1 ≤ k ≤ s such that j ∈ G k . Then, by (S.10), , which would imply by (S.12) that Ψ λ (x, u) = x + v r for some r ∈ G k,− i . Since we are assuming that Ψ λ (x, u) = x + v ℓ and we know the vectors v 1 , . . . , v n are distinct, we obtain that ℓ = r ∈ G k,− i . This contradicts ℓ / ∈ G k,− i .
ii) if ⟨A i• , v j ⟩ = 0 and ⟨A i• , v ℓ ⟩ = 1, then j / ∈ G k,+ i and ℓ ∈ G k,+ i . By (S.7), we would then have r∈G k,+ iΥ r (y) ≥ r∈G k,+ i Υ r (x), which would imply by (S.11) that j ∈ G k,+ i . This contradicts j / ∈ G k,+ i . iii) in all the other cases, that is when Combining the above for case a), we see that To see this, observe that for every i / ∈ K y , ⟨A i• , y − x⟩ > 0 and as for (S.16), which would imply by (S.12) that Ψ λ (x, u) = x + v ℓ for some ℓ ∈ G k,− i , but this contradicts the assumption that Ψ λ (x, u) = x. So we must have Case 2:Ψ λ (y, u) = y. Again, we consider two subcases. a) If Ψ λ (x, u) = x, then (S.8) holds, because x ≼ A y. b) If Ψ λ (x, u) = x + v j for some 1 ≤ j ≤ n, we claim that y ∈ K A + x + v j for the corresponding value of j. To see this, fix the value of j for which Ψ λ (x, u) = x + v j , let 1 ≤ k ≤ s be such that j ∈ G k , and observe that for every i / ∈ K y , ⟨A i• , y − x⟩ > 0 and as for (S.17), i , then by (S.7), we would have ℓ∈G k,+ iΥ ℓ (y) ≥ ℓ∈G k,+ i Υ ℓ (x), which would imply by (S.11) thatΨ λ (y, u) = y + v ℓ for some ℓ ∈ G k,+ i , but this contradicts the assumption thatΨ λ (y, u) = y. So we must have ⟨A i• , v j ⟩ ≤ 0 and hence ⟨A i• , y − (x + v j )⟩ ≥ 0 for all i ∈ K y . Thus, we have y ∈ K A +x+v j and then Ψ λ (x, u) = x+v j ≼ A y =Ψ λ (y, u).
In order to prove Theorem S.2, from here on we can follow a similar procedure to the one used in the proof of Theorem 3.1 after Lemma 5.1 was proved there. For the case where (5.1) holds, we define two discrete-time processes, Y = (Y k ) k≥0 andY = (Y k ) k≥0 , by defining Y 0 := x • ,Y 0 :=x • , and for k ≥ 0, (S.20) and define X andX using these and an independent Poisson process N as in (5.13). For the case where (5.1) does not hold, we can use a truncation procedure similar to that for Theorem 3.1. In both cases, we use Lemma S.1 instead of Lemma 5.1.